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Abstract 

The double-layer Heisenberg antiferromagnet with intra- and inter-layer cou- 
plings J and Jj_ exhibits a zero temperature quantum phase transition be- 
tween a quantum disordered dimer phase for g > g c and a Neel phase with long 
range antiferromagnetic order for g < g c , where g = Jj_/J and g c ~ 2.5. We 
consider the behavior of the system at finite temperature for g > g c using two 
different and complementary approaches; an analytical Brueckner approxima- 
tion and numerically exact quantum Monte Carlo simulations. We calculate 
the temperature dependent spin excitation spectrum (including the triplet 
gap), dynamic and static structure factors, the specific heat, and the uniform 
magnetic susceptibility. The agreement between the analytical and numerical 
approaches is excellent. For T — > and g — > g c , our analytical results for the 
specific heat and the magnetic susceptibility coincide with those previously 
obtained within the nonlinear a model approach for N —* oo. Our quan- 
tum Monte Carlo simulations extend to significantly lower temperatures than 
previously, allowing us to obtain accurate results for the asymptotic quan- 
tum critical behavior. We also obtain an improved estimate for the critical 
coupling: g c = 2.525 ± 0.002. 

PACS numbers: 75.10.Jm, 75.30.Kz, 75.50Ee 



Typeset using REVTpjX 



I. INTRODUCTION 



The Heisenberg antiferromagnet on a double-layer square lattice has become an impor- 
tant model for studying quantum antiferromagnetism in two dimensions. It was introduced 
by Millis and Monien as a phenomenological model to capture the spin-gap behavior ob- 
served in high-T c compounds such as YBa2Cu3O6+ x ,0B in which the basic structural unit is 
a pair of CuC>2 planes. The Heisenberg bilayer with inter- and intra-layer couplings J and 
J± can be tuned through an antiferromagnetic order- disorder transition at zero temperature 
by varying the ratio g = J±/J. For g < g c there is long-range antiferromagnetic order at 
zero temperature, whereas for g > g c the tendency to singlet fomation across the layers leads 
to a disordered ground state and a spin gap. The critical ratio g c ~ 2.5 has been deter- 



mined using a variety of numerical methods.ETH The inter-layer coupling in bilayer cuprates 
is typically much smaller than this value; neutron scattering experiments indicate g ~ 0.1 in 
YBa 2 Cu 3 06S Although experimental realizations of antiferromagnetic bilayers with g ~ g c 
are currently lacking, the model in this regime is important for theoretical investigations of 
quantum critical and quantum disordered behavior.0 

It has been argued by HaldanJ and Chakravarty, Halperin and Nelson! that the T = 
quantum phase transition in two-dimensional quantum antiferromagnets is described by the 
2+1-dimensional nonlinear 0(3) <r-model, and therefore the universality class should be that 
of the 3D classical Heisenberg model. The nonlinear a-model has three distinct regimes in 
the T — g plane near the quantum critical point. For g < g c the long-range antiferromagnetic 
order at T = is destroyed by thermal fluctuations for any T > 0. At low temperatures, in 
the so-called renormalized classical regime, the correlation length diverges exponentially as 
T — ► 0; £ oc e 2?rps//T , where p s is the spin stiffness. For g > g c the singlet-triplet excitation 
gap implies a low-temperature quantum disordered regime where the correlation length is 
temperature independent and remains finite as T — > 0. Exactly at g c , p s vanishes and the 
correlation length diverges as 1/T. This is the leading behavior also in the high-temperature 
quantum critical regime for both g < g c and g > g c . Exactly at g = g c , the quantum critical 
regime extends down to T = 0, whereas for g ^ g c there is a crossover to either renormalized 
classical or quantum disordered behavior at lower temperature. 

The mapping of the Heisenberg model to the nonlinear a-model has been rigor- 
ously proven only for the single-layer square lattice antiferromagnet with nearest-neighbor 
interactions,! which itself does not exhibit a zero-temperature quantum phase transition. 
Assuming the cx-model description, the universal dynamic and static properties of two- 
dimensional antiferromagnets in the vicinity of a zero-temperature phase transition have 
been studied in detail by Sachdev, Ye, and Chubukov0'0 using an 1/iV expansion (N = 3 
is the number of components of the a-model). They found that close to criticality, many 
physical observables such as the specific heat, the magnetic susceptibility, e.t.c, depend in 
a universal manner on a small number of model dependent parameters. The Heisenberg 
bilayer is an ideal model to numerically test these predictions.u The model is not frustrated, 
and therefore the "sign problem" hampering simulations of the single-layer J1-J2 model is 
not present. Compared to other dimerized models, the bilayer has the advantage that the 
square lattice symmetry in the planes is not broken, and therefore the asymptotic behavior 
should be more easily observed.il 

The zero-temperature transition of the bilayer system has been well studied by various 
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numerical methods, including dimer-series expansions by Hid J and Zhengjj and quantum 
Monte Carlo simulations by Sandvik and Scalapino.! These studies determined the critical 
point in the range be g c m 2.50 — 2.56. The scaling behavior is consistent with 3D Heisenberg 
critical exponents. 

The finite temperature properties have been also extensively studied by numerical meth- 
ods. Quantum Monte Carlo simulations were performed by Sandvik and coUaboratorsilli 
and showed a reasonably good agreement with the predicted quantum critical behavior. 
High-temperature expansions were carried out by Oitmaa et alM and finite-cluster exact 
diagonalisation was used by Jaklic and Prelovsek.EJ Finally, high-order strong-coupling ther- 
modynamic expansions for this model were developed by Elstner and Singh.li3 Unfortunately, 
all the methods have so far been restricted to relatively high temperatures. A particular 
recent concern is that the strong-coupling expansions by Elstner and Singh indicate a change 
in the behavior of the susceptibility and the specific heat close to the lowest temperature 
considered in the Monte Carlo simulations.^ It is not clear whether this is due to the break- 
down of the combined high-temperature and strong-coupling expansion, which is expected 
at low temperature, or if the Monte Carlo results do not reflect the true asymptotic quantum 
critical regime. 

On the analytical side, until recently the progress with the bilayer model was less im- 
pressive because the conventional spin-wave theory used by HidJ and the Schwinger boson 
mean field approach applied by Millis and Monien0 were not able to take into account in 
an appropriate way the strong interactions between elementary triplets. These approaches 
yielded results which were inconsistent with the numerical studies. In particular, a much too 
high critical coupling g c m 4.5 was obtained. A more sophisticated treatment by Chubukov 
and Morr,0 taking into account longitudinal spin fluctuations which are important close 
to g c , gave an improved value g c « 2.7 and was also able to explain other aspects of the 
quantum Monte Carlo results. 

A very efficient analytical approach to deal with a dimerized disordered quantum spin 
system at zero temperature has recently been developed by Kotov et al. .El This method is 
based on the same ideas as the Brueckner theory of nuclear matter,@ and we therefore call 
it the "Brueckner method" . This approach has already led to rather spectacular success in 
several problems, including a calculation of g c (yelding g c = 2.60}. the zero temperature crit- 
ical index, and the excitation spectrum for the bilayer model,0@ the elementary excitation 
spectrum and bound states for one-dimensional spin-ladders and spin chains§3'i3, as well as 
the ground state structure and the excitation spectrum of the frustrated two-dimensional 
J i — J 2 model.il 

One of the purposes of the present work is to develop a generalization of the Brueckner 
method to the case of finite temperature, using as a testing ground the Heisenberg bilayer 
Hamiltonian. All analytical results are general, however, and can be applied to any dimerized 
quantum spin system, including one-dimensional spin ladders and chains. In the present 
paper we focus on the quantum disordered and quantum critical regimes of the bilayer. Our 
formulas can be formally extended to the renormalized classical regime, but in that case 
there are very important corrections that we plan to consider in a future study. 

In order to test the reliability of this new theoretical approach, we have also carried out 
quantum Monte Carlo simulations at considerably lower temperatures than before. This al- 
lows us to definitely conclude that the asymptotic quantum critical regime has been reached, 
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and to resolve the above mentioned questions arising from strong-coupling expansions. The 
simulations also give an improved estimate of the critical coupling; g c = 2.525 ± 0.002. 

Using the Brueckner method, we have calculated the temperature dependent spectrum 
of excitations, dynamic and static structure factors and susceptibilities, the specific heat, 
and the uniform magnetic susceptibility. In the limit T — » 0, g — > g c , our results for the 
correlation length, the specific heat, and the magnetic susceptibility coincide with those 
previously obtained using N — > oo calculations for the nonlinear a-mo del.0'0 The Monte 
Carlo approach gives numerically exact results for the thermodynamics. The agreement 
with the analytical calculations is excellent at low temperature (T/J < 0.5). Using sum 
rules, we can also extract an approximate excitation spectrum from the quantum Monte 
Carlo results. Also in this case the agreement with the theory is very good. 

The rest of the paper is organized as follows. In Sec. II we first summarize the main steps 
of the Brueckner approach at zero temperature, and then proceed to the new calculations 
at finite temperature. In Sec. Ill we discuss the Monte Carlo simulations and present data 
allowing us to extract an accurate value for g c . We also discuss finite-size effects. We 
present comparisons between analytical and numerical results for various physical quantities 
in Sec. IV. In Sec. V we summarize and discuss our main conclusions. 



II. BRUECKNER THEORY 

The double-layer Heisenberg model we study is defined by the Hamiltonian 

H = J 2j(Sii ■ Sij + S 2 j • S 2 j) + J± Sii ■ S 2 j, (1) 

where S a i is a spin-1/2 operator at site % of layer a (a = 1, 2), denotes a pair of nearest- 
neighbor sites on the square lattice, and the total number of sites in a layer is L 2 . Both 
the coupling constants are antiferromagnetic, i.e., J, J± > 0. In Sec. II-A we summarize the 



main steps of the Brueckner diagrammatic approach developed in Refs. p!8|j20| for this model 
at zero temperature. In Sec. II-B we generalize the approach to T > 0. 



A. Zero temperature 

Our approach is formulated in the basis of bond singlets and triplets. We define a creation 
operator t* ai for a triplet (S = 1) with polarization a — x, y,z at bond i, where % connects 
two nearest-neighbor spins Si(z) and S 2 (2). In the bond operator representational! 

Sx,S) = \{±t a ±A- ie^ty), (2) 
and one can exactly map the Hamiltonian ([!]) to the effective Hamiltonian 

#eff = H 2 + H 4 + H U} (3) 

where 
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H 2 — ^ J±taitai + ~^ (*L*«+1 + *L*L+1 +h.c), (4a) 



a,, 
J 



Ha — (tcdtpjtpitaj - taitltjtpitpj), (4b) 

(i,j),a,P 

Hu = UY tlitpM^ U 00 ' ( 4c ) 

Hereafter we set A = 1, except at the end of Sec. II-B where it will be convenient to use a 
variable A to formally separate the contributions from quadratic and quartic terms in the 
Hamiltonian. 

At the quadratic level the Hamiltonian @ can be diagonalized by a combination of 
Fourier, 

C = ^E4, Q ^ k+k ^ k = (-,-), (5) 

and Bogoliubov, 

tka = u k a ka + v k a ] _ ka , (6) 
transformations. This gives the excitation spectrum 



2 -K- Bl (7) 



k 



where 



A k = J ± + 2 Ja, (8a) 

B k = 2J£ k , (8b) 

and 

£k = — (cosfc x . + cosk y )/2. (9) 

The momentum takes values within the Brillouin zone — ir < k x , k y < n, but for convenience 
we have shifted the argument in the Fourier transform (^) by ko = (ir, 7r). In this notation 
the minimum of the spin- wave dispersion is at k = 0. 

An infinite on-site repulsion between triplets, Hu, has been introduced in order to take 
into account the hard-core constraint = (only one triplet can be excited on a bond). 
For g > g c m 2.5, the contribution of the quartic term if 4 is relatively small and Hy therefore 
gives the dominant contribution to the renormalization of the spin-wave spectrum. Since 
Hjj is infinite the exact scattering amplitude 



a/3, 7<5 



T(K)(5 ai 5 p& + 5 a& 5 M ) } (10) 



where K = (k, uj) is the total energy and momentum of the incoming particles, can be found 
from the following equation, which is diagrammatically represented in Figure [11(a), 

r(K) = -(4£ u2pULp ) 1 ■ (ii) 
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The Bogoliubov coefficients are 

ulvl = ±l/2 + A k /(2Lu k ). (12) 

The basic approximation made in the derivation of T(K) is the neglect of all anomalous 
scattering vertices which are present in the theory due to existence of the anomalous Green's 
function, 

G a (k,t) = -^(T(t t _ kQ (t)tL(0))). (13) 

The crucial observation is that all anomalous contributions are suppressed by a small pa- 
rameter of the theory — the density of triplet excitations, 

p = (iu = ^E< (14) 

It was found in Ref. [18| that p ~ 0.12 at the critical point, g = g c , and that p decreases as g 
increases. 

The normal self-energy corresponding to the scattering amplitude T(K) is given by the 
equation shown diagrammatically in Figure 0(b): 

E(k, U ) = -i J2 <r(k + q, u - Wq ). (15) 

In order to find the renormalized spectrum, one has to solve the coupled Dyson equations 
for the normal Green's function 

G n (k,t) = -i(T(t ka (t)4 a (0))) (16) 

and the anomalous (|T3|). Further separation into a quasiparticle contribution and incoherent 
background yields 



where 



u k = Z^Al-Bl (17) 



A k = J ± + 2J& + S(k, 0) + 4 ]T £ q <, (18a) 

L q 

5 k = 2 J£ k - 4J£ k -i £ f qUqUq . (18b) 



We also define 



<9E x 1 



(19) 



and 



Ui,Vi = ^t4 (20) 
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The renormalized coefficients and fll8|) also take into account the quartic interaction 
( fi"b"D in the one-loop approximation because the corresponding effect is small. In order to 
find the spectrum, equations (|TTD , (|i5|) , and ([L7])-(^(]) with the substitutions 



«k -> V^k^k, ^k -> V^fcVk (21) 

have to be solved self-consistently for £(k, 0) and Z^. 

The results of the self-consistent solution for the spin-wave gap, A = w&=o, are in 
excellent agreement with a dimer series expansions^ even close to the critical point, as 
shown in Figure 0. The value of the critical coupling is g c ~ 2.60, only slightly larger than 
the numerically determined value g c = 2.51 ± 0.02.Q It has been demonstrated in Ref. ^0 
that the Brueckner approximation itself results in a critical behavior for the triplet gap, 
Ao ~ {g — 9c) u , with critical index v — 1/2. It was also shown that for the model under 
consideration the actual small parameter is pln(J/A ). Thus when the gap is very small the 
approach can fail. Nevertheless, the critical behavior of the gap can be analyzed analytically 
(see Ref. ^) by expanding the Brueckner equations in powers of the triplet density p. In the 
leading approximation the contribution from the quasiparticle residue to the renormalization 
(PT|) can be neglected and we have to set u\ — 1 in the vertex (|TT|): 

r(q, - Wq ) = ( I ( d X - \ ) 1 . (22) 

Then the variation of the self-energy flTB]) with the deviation 5g = g — g c results in a linear 
relation between gap, A , and 5g near the critical point: 

A /Jna(g-g c ), (23) 

where the constant a ~ 1.1 has been evaluated in Ref. [2(]. The critical exponent v = 1 
corresponding to (|23D is the same as the one predicted within the nonlinear a model for 
iV — f oo.0 The plot of the gap, Ao(g), obtained from the self-consistent solution of low- 
density Bruekner equations, i.e., the solution of Eqs. (p~5[ ) and (p^)-([20|) together with the 



vertex ( p2|) [and the substitutions u^, v k — > U^, Vu\ is also shown in Figure |. The resulting 
value of the critical coupling, g c w 2.32, is about 10% lower than numerically determined 
one. 

Consideration of the first correction due to the finite triplet density, i.e. using uiv? 



p u k-p 



v. 



p)(l + ^k-p) ~ 1 + + fk_p in the vertex (0), results in a logarithmic correction at 



small momenta, T(q, —oo q ) ~ T c (l + /Sing), and yields 

Ao/J = a6g(l-phi6g), (24) 



with the constant (3 ~ 0.3 evaluated in Ref. ^0| Assuming scaling behavior, Aq ~ {&9) v ', the 



logarithmic correction results in a critical index, v = 1 — /3 « 0.67 ±0.03, which agrees with 
numerical studies and the prediction of the a model [y ps 0.70). This concludes the solution 
of the two-layer Heisenberg antiferromagnet at zero temperature. 
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B. Finite temperature 



Having described in the previous section the Brueckner method at zero temperature as a 
starting point, we next consider the bilayer Heisenberg antiferromagnet at finite temperature. 
The generalization of the method to T > is simple because in essence it is a low-density 
approximation. Thus it is convenient to use the Feynman technique instead of the Matsubara 
one. In this approximation they are equivalent to each other. 

Only the normal Green's function G n (k,t) = — 2(T(a kQ (i)a ka (0))) exists for physical 
operators a kce . At finite temperature in an idial gas approximation, it can be written via 
retarded and advanced functions, G R,A (k,u) = [uj — uj k ± «0] _1 , as (see Ref. |2B|) 

G B (k,w) = \{G R + G A ] + ictanh^G* - G A \. (25) 

Then, after Fourier @ and Bogoliubov @ transformations the normal and anomalous 
Green's functions 

Gl(k,t) = -i{T(t ka (t)tl a (0))), (26a) 
G2"(k,t) = -i<r(ti te (t)tL(0))>, (26b) 

at finite temperature are 

g t = uj(l + n k ) _ uln k _ vl(l + ra k ) + vjn^ ^ 

11 ' uj — cj k + iO uj — co> k — iO uj + aj k — iO uj + a; k + iO ' 

wr/i \ / l+n k n k n k 1 + n k \ 

G a (k, u) = u k v k — — + — — — , (28) 

\uj — u; k + i\} uj — LUk — tO uj + + tO uj + uJk — tOJ 

where n k = (a kQ ak,a) = (e^ k / T — 1) _1 with no summation over a. The essential feature is 
that despite the bosonic commutation relations of the single triplet operators, these do not 
have conventional bosonic distribution function because of the strong interactions. Actually, 
aj k is a functional of n k , and thus the expression n k = (e Wk / T — 1) _1 should be considered 
as an implicit definition of the thermal distribution function n k . The density of triplet 
excitations, trivially obtained from the normal Green's function ( p7|) at finite temperature, 
is given by 



P 



T ~ (40 = i ■ r lim G T n (k, ")e— d £ 

4W + n J+^ ( 29 ) 

Our approach is valid for low density of triplet excitations, and therefore our consideration is 
restricted to low temperatures; T < J±, J. We will take into account only leading tempera- 
ture corrections to the self energy (|T5|) neglecting temperature corrections to the vertex (p2|). 
The quantum transition disappears at any finite temperature and there is only a disordered 
phase with a finite gap towards to triplet excitations. Thus we expect that the low-density 
Brueckner equations is an appropriate tool for the model in the g > g c region, where the 
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neglected logarithmic corrections, pin (Ay /J) and pln(Ao/J) are not large and the vertex in 
the form ( p2|) can be used. Using the normal time-dependant Green's function ( P7|) for the 
diagram yields 

£ T (k, u) = + ^q)r(k + q, u - c^ q ) + <n q r(k + q, u + w q )l . (30) 

In order to find the renormalized spectrum at finite temperature, one has to solve the coupled 
Dyson equations for the normal and anomalous Green's functions (p7[) and (p8[) , as was done 
in Sec. II-A at zero temperature. Further separation into a quasiparticle contribution and 
incoherent background renormalizes the coefficients 

A k = Jl + 2J£ k + S T (k, 0) + 4J£ k -^ X! C q [<(l + n q ) + <n q ], (31a) 

^ q 

5 k = 2J£ k - 4J£ k -^ J] &Wq(! + 2n q ), (31b) 

^ q 

where the effect due to H4 is also taken into account in the the one-loop approximation. Then 
the finite temperature renormalized spectrum u; k , the renormalization constant Z k , and the 
Bogoliubov coefficients ?7 k and Vk are given by Eqs. flTT] ) and ( PUD with the self-energy ( |3"0"| ) 
and the coefficients ([Hp. 

In the low-density approximation the equations (p!7|),(p!9D,(|20D,(p2]),(pO|), and (|3T| ) should 
be solved self-consistently for S T (k, 0) and Z(k) with the substitution w k — > f/ k , f k — > Vk- 

All thermodynamic functions of the system can be obtained from free energy F. The 
total energy of the system can not be found as a simple summation of quasiparticle energies 
because of the strong interactions, but it can be calculated in following way: Assume that 
the parameter A in the quadratic term ([|a|) takes values between and 1. Then we have 



dF _ /dH\ 
~d\ ~ \~d\ / T ' 



(32) 



If A = 0, the system is a set of dimers with free energy 

F = -Tln(l + 3e- J±/T ). (33) 
Integrating Eq. (|32| ) then yields 



where 



/ dH \ J 1—-^ j- j| 

( ~qT~ ) = ~2lJ2 £^,\t a it a i+l + t a it ai+ i + h.c.) 

6 

= J TI E + nj + u 2 k n k + u k v k [l + 2n k ]). (35) 

^ k 

Using Eq. (]34]), one can calculate the thermodynamic potential and all thermodynamic func- 
tions using the self-consistent solution for the spectrum tu k and the Bogoliubov coefficients 
t>k and w k at finite temperature. 
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This concludes the solution of the two-layer Heisenberg antiferromagnet model at finite 
temperature. The results of self-consistent solution of the Brueckner equations for the triplet 
gap, the spectrum, as well as various thermodynamic quantities will be presented in Sec. IV, 
along with comparisons with numerically exact results obtained using the quantum Monte 
Carlo method discussed in the next section. 



III. QUANTUM MONTE CARLO 

Quantum Monte Carlo simulations of the Heisenberg bilayer have previously been used 
to determine the T = critical coupling, with the result g = 2.51 ± 0.02.Q This value 
is also in good agreement with results obtained using dimer series expansion methods.!'! 
Finite-temperature simulations in the quantum critical regimeSE! have shown generally good 
agreement with predictions of 1/N calculations for the nonlinear a model. However, the al- 
gorithm used in the previous simulations did not work well at very low temperatures, and 
therefore there has been some concernEl that the results do not reflect the true T — > be- 
havior. Recently, the "stochastic series expansion" algorithm0@ used in previous work has 
been considerably improved by the introduction^ of a new cluster-type updating scheme 
(inspired by the "loop algorithms" developed for worldline Monte Carlo simulations!^) that 
overcomes the limitations of the previous local sampling scheme. We have used this algo- 
rithm for large-scale calculations in order to obtain a more precise estimate of the critical 
point and to study the low-temperature quantum critical and quantum disordered behavior. 
In this section we extract the T = critical coupling and also discuss effects of finite lattice 
size in T > calculations. In the Sec. IV we will compare results for the thermodynam- 
ics and the excitation spectrum with the Brueckner theory. For details of the simulation 
algorithm we refer to previous literature.il'i! 

One way to extract the critical point is from the size dependence of the T = spin 
stiffness p s . Imposing a uniform twist such that the spin-spin interaction is modified 
according to 

SZ£% + S y ai Sl 3 -> + SISD cos (0) + {SlS* aj - S^) sin (0) (36) 

(for all nearest-neighbor spin pairs in either the x or y lattice direction in both layers a = 1, 2) 
the stiffness is given by 

Ps = (37) 

where E((p) is the internal energy per spin. The stiffness can be related to the fluctuation of 
the "winding number" in the Monte Carlo simulationsil and can hence be evaluated directly 
without actually including the twist. 

We study quadratic lattices with linear size L, i.e., the total number of spins is 2L 2 . 
According to scaling theory,!! at the critical point p s should depend on the system size 
according to 

p s (q = g c ) ~ L d ~ 2 ~\ (38) 
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where d is the dimensionality and z is the dynamic exponent. Hence, with z = 1, as follows if 
the transition is in the universality class of the 3D classical Heisenberg model, Lp s graphed 
versus g for different system sizes should intersect at g c . In Figure ^ we show results for 
L up to 20, calculated at temperatures sufficiently low to give the ground state (for the 
largest system, L = 20, J/T = 320 was used). For the largest lattices, the curves intersect 
at g ~ 2.525, in good agreement with the previous Monte Carlo estimate g c = 2.51 ± 0.020 
Other useful quantities that can be calculated in the simulations include the static struc- 
ture factors 5 ,± (q) and susceptibilities X (q)) which are defined according to 



S± (^ = 4l E e 4q - (r2 - ri)) [Si, n ± S 2>ri ][S ltr2 ± S 2 , r2 ] (39a) 

lL ri,r 2 

X ± (q) = E e iq - (r2 - ri)) [dr[S ltri {r) ± S 2 , n (r)] (0) ± S 2 , r2 (0)]. (39b) 



ri,r2 



We first consider the uniform magnetic susceptibility, x — X + (0, 0), close to the estimated 
g c . According to calculations for the nonlinear a model, x{T) should be linear in T in the 
quantum critical regime, with intercept at g = g c . In order to study the behavior at low 
temperatures, lattices sufficiently large to eliminate finite-size effects have to be used. In 
Figure [| we graph results at g = 2.53 for several system sizes. The results show that in 
order to obtain results reflecting the thermodynamic limit down to T/J = 0.1, systems with 
L = 64 — 128 are required. Note that the asymptotic T — > behavior for a finite lattice 
is always an exponential decay to zero, reflecting the finite-size singlet-triplet gap (which 
scales to zero with increasing L for g < g c and remains finite for g > g c ). However, in the 
temperature regime shown in Fig. |4], the finite-size effects actually enhance x- 

In Figure |5] we show results for system sizes sufficiently large to eliminate finite-size 
effects for g = 2.52 and g = 2.53. For both couplings, the behavior is very close to linear for 
0.1 < T/J < 0.4. The intercept of a line fitted to L = 128 data for T/J < 0.17 is positive 
for g = 2.52 and negative for g = 2.53, in consistency with g c w 2.525 estimated from the 
T = stiffness scaling. Using the data shown in Fig. we estimate g c ~ 2.525 ±0.002. This 
result is slightly lower than estimates (g c = 2.54 — 2.56) from dimer series expansions.!'! 

In Sec. IV we will make use of sum rules and the ratio 

^w = a (40) 

to extract the approximate finite-temperature triplet excitation spectrum. Here, in Figure || 
we show the finite-size effects in R~(ir, 7r) for g = 2.52 and 2.53. The classical value of 
i? ± (q) = 1 for any q, and this is also expected to hold for R~(tt,7t) of an antiferromagnet 
in the renormalized classical regime (where both the staggered structure factor and the 
staggered susceptibility diverge). For a system with a gap, on the other hand, both 5' ± (q) 
and x ± (q) converge to finite values for all q and R~(TT,n) therefore diverges as 1/T at low 
temperature. Nonlinear cx-model calculations have predicted a temperature independent 
R~(tt,ti) in the quantum critical regime, with the value i?~(7r,7r) w 1.09.00 Previous 
Monte Carlo calculations showed a reasonable agreement with this prediction close to the 
critical point for T '/ J m 0.3.0 The results shown in Fig. [] demonstrate that the ratio is very 
sensitive to the system size at lower temperatures. The finite singlet-triplet gap present 
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for any finite L causes a 1/T behavior of R as T — > 0, starting at a temperature that 
decreases with increasing L. The results shown in Fig. ||| indicate that L = 128 should be 
sufficiently large to eliminate finite-size effects for all T/J > 0.1. For this lattice size, the 
data for g = 2.52 decreases with decreasing T whereas there is an increase for g = 2.53. 
This behavior is consistent with the predicted temperature independence exactly at g c and 
g c 2.525. In Sec. IV we will present analytical results for the temperature corrections to 
i?~(7r,7r) and make further comparisons with numerical results. 



IV. RESULTS 

We now discuss finite-temperature calculations at the critical point and in the quantum 
disordered regime. Results of the low-density Brueckner theory will be compared with 
quantum Monte Carlo in the g > g c region. Note that the critical point resulting from 
the self-consistent solution of the low-density Brueckner equations at zero temperature, 
g c w 2.320 (see Figure 2), is different from the actual value g c ~ 2.525. Unless stated 
otherwise, we will consider numerical results at g c as the average of data for g = 2.52 and 
2.53, which should be close to the actual g c ~ 2.525 critical behavior in the temperature 
regime considered. We compare the results with the Brueckner theory predictions for the 
critical point in that approach; g = 2.320. To compare our analytical and numerical results in 
the quantum disordered regime we will take g = 3.0 (where the Brueckner theory predictions 
are in perfect agreement with the numerical data at zero temperature) both for the Brueckner 
and Monte Carlo calculations. 

A. Triplet gap 

Let us first consider the critical behavior of the triplet gap, A^ = Uk=o, described by 
the Brueckner equations. The temperature and the deviation from the critical point must 
satisfy the condition T / J, Sg <C g c . No restriction is imposed on the ratio Sg/T, however. 

Close to the critical point and for small momenta (fcc 1) the dispersion can be repre- 
sented as 

u k « ^A* + 7 2 P, (41) 

where 7 = 1.9J is the zero temperature spin-wave velocity.iSil Eq. (|TT| ) at the point k = 
gives 

A 2 T = Z (Al-B 2 ). (42) 

It is convenient to introduce the values of A ,B , and Z at zero temperature 
A T=0 , B T=0 , Z T=0 , which satisfy Aq = Z T=0 (A^ =0 — B^ =0 ). Let us vary the temperature, 
keeping J± and J fixed. The main contribution to the variations of the integrals in Eqs. ( |3Tj| ) 
and ( |3"TD comes from small momenta; g ~ A/7 < 1. Then we find the variations of B k 
at k = to leading order: 
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Ao = A T =o + 6Z T '(0,0) - J -^R{T) 



B = B T=0 + 



JA C Z C 

7T7 2 



R(T), 



where 



5s T (o,o)«^i?(r), 



and 



J2(T) = -A, + Ao + 2T /~ 



(43a) 
(43b) 

(44) 

(45) 



and where r c , Z c and A c denote the vertex T(0,0), the quasiparticle residue Z k=0 , and A 
evaluated at T = and g = g c . 

If we substitute ( fi3"D into ( f4*2"|) and neglect terms quadratic in At, A , we find that the 
variation of the self-energy has to vanish; <5£ T (0,0) = 0. Thus R(T) = determines the 
finite-temperature gap A T as a function of the zero temperature gap A : 



Also, we have found that 



At _ Ao f 00 dx 
T ~ T J^z. e x - 



5p T « ^i?(T) = 0. 



(46) 



47T7 Z 



(47) 



Thus, as the temperature increases the quantum fluctuations are reduced while the temper- 
ature fluctuations are enhanced, the total triplet density remaining constant. After a simple 
integration of fl46|) we obtain the analytical solution 



At = Ty(x) 



where 



y(x) 



21n 




= 2arcsh 






.2 



(48) 



(49) 



with x = Ao/T. The variable x determines whether the antiferro magnet is in the quantum 
critical (x — *■ 0), quantum disordered (x — > oo), or the renormalized classical (x — ► — oo) 
region. It corresponds to the variables x\ = —2 / x and x^ — l/x introduced within the non- 
linear a model formalismS for the critical regions in the T — g plane. The zero temperature 
gap, A / J = a(g — g c ), is defined for positive values (g > g c ) but we analytically continue it 
to negative values (g < g c ). The temperature gap A T depends on the coupling g via the zero 
temperature gap A (g). Thus, in order to calculate A T (g) one should substitute Eq. 
into (H). 
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The universal function y(x) for the inverse correlation length, l/£ = At, is exactly the 
same as that previously obtained for the nonlinear a model (N — > oo),0 with the limiting 
behavior 

, 1 + V§ x 

y{x) = 21n^^ + -j=, x^O, (50a) 

y(x) « x + 2e _:r , x — ► oo, (50b) 
y(x) « e^ 2 , x — > -oo. (50c) 

Thus the correlation length is exponentially large in the renormalized classical region, pro- 
portional to 1/T in the quantum-critical region, and of the order 1/Aq in the quantum 
disordered regime. At the crititical point (x = 0) we have 



A T = yoT, y = 21n 



1 + V5 



0.962424. (51) 



In Figure [7| we compare the analytical result (f48|) at g — g c , and g = 3.0 obtained under the 
assumption of small T and Sg with the corresponding self-consistent numerical solutions of 
the Brueckner equations for the triplet gap. The agreement is very good for g = g c at low 
temperatures, and slightly worse at g — 3.0. 

Now, let us consider the following dimensionless ratio 

*(q) = P^-y (52) 

^x(q) 

where S(q), are the odd (— ) static structure factor and susceptibility defined in 

Eqs. (|39| ) [with the momentum shifted by (ir, 7r)]. They are related to the corresponding 
dynamic structure factor, 

S(q,w) = Im| e i " t <T(S(q,*)S(-q,0)))dt, (53) 

where Sj = — S 2 j, via the sum rulesH 

S(q) = - / SYq, w) ( 1 + e^ /T J dw, (54a) 

7T JO V 7 

X (q) = ?f % U )(l-e-^. (54b) 

7T JO V 7 W 



Taking into account the contribution from the elementary triplet only, the structure factor 
can be written in the form 

S(q,u) = (w q + w q ) 2 <5(^-u; q ). (55) 

Then, using the above sum rules ( p4a| ) and ( |54b| ), the dimensionless ratio -R(q) can be 
rewritten in the form 

1 -|_ p-^q/ r ij 
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The importance of this ratio is that it can be obtained using Monte Carlo simulations 
by calculating the static susceptibility and structure factor, as discussed in Sec. III. The 
spectrum of elementary excitations, c<j q , can then be extracted from Eq. (|5(J). At zero 
momentum, q = 0, the critical behavior of R{0) is given by 

Then the limiting behavior of -R(O) is readily obtained from Eq. fll^p: 

fl(0) = ^(yo + ^(l- A=Vo)), (58a) 
R(0) = | + e - x (l + x), x^oo, (58b) 

x/2 

R(0) = 1 - x -> -oo. (58c) 

We have calculated R(0) as a function of temperature using quantum Monte Carlo sim- 
ulations close to the critical point, at g — 2.52 and 2.53, and using the Brueckner equations 
at the critical point (x = 0). As discussed in Sec. Ill (Fig. f|), R{0) is very sensitive to 
the lattice size at low temperatures. In Figure |] we show Monte Carlo results for systems 
sufficiently large to eliminate finite-size effects (linear size L up to 128). In Sec. Ill the 
critical point was determined to be g c = 2.525 ± 0.002. Even for \g — g c \ as small as 0.005, 
there are very large differences between the g — g c > and g — g c < data. In order to 
extract the behavior exactly at g — g c , simulations would have to be carried out for values 
of g also between 2.52 and 2.53, and the critical point would have to be determined to even 
higher accuracy. The available data nevertheless shows that the Brueckner theory is very 
accurate. The analytical curve falls between the two numerical data sets at low temperature, 
and the T dependence at high temperatures is also very well reproduced. In the Brueckner 
theory R{0) — > 1.076 at g — > g c as T — > 0. The numerical results suggest that the exact 
value is slightly higher, but close to 1.08. 

The temperature dependent triplet gap extracted from the Monte Carlo data for R(0) 
using Eq. ( |56| ) is presented in Figure [7|, along with our analytical prediction. The agreement 
is almost perfect at the critical point, while at g — 3.0 it is not so good but still quite 
reasonable. 

The formula (|56|) is not exact, due to the neglected contribution of the three-particle 
continuum and broadening. Note that the contribution of the two-particle continuum is 
canceled by symmetry. The three-particle continuum is negligible only near the q = (0, 0) 
at the critical point while it can be important in other cases. 

Using Eq. fl5"6"D , we can extract the triplet spectrum u>k from the Monte Carlo results for 
.R(q). In Figure |9| we compare it with the solution of the Brueckner equations at g = g c and 
T = 0.5J. The agreement between the Monte Carlo and analytical results is even better 
then might expected (the agreement remains very good also at lower temperatures). The 
Brueckner approximation should be accurate if the triplet density p is low, but in fact it is 
as high as p ~ 0.1.0 The surprisingly good agreement justifies Eq. (|56D and the assumption 
that the dynamic structure factor is peaked at u = oJk [see Eq. flSSp]. Thus the broadening is 
negligible and the elementary triplet is well defined for g > g c at low temperature. However, 
the broadening can be important for g < g c . 
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B. Specific heat 



The calculation of the specific heat, Cy, is straightforward. The energy of the system 
cannot be written as a simple summation over quasiparticle energies, E ^ Y^^k^k, due to 
the strong interaction. However, one can write the energy variation as 5E = Y^k^k, i n 
analogy with the Landau theory of Fermi liquids. Thus, the specific heat per spin is 

* = 2Z*?<*fr (59) 

k 



It is convenient to use Eq. (|59|) instead of taking of the second derivative, Cy = 
— ^Td 2 F/dT 2 , using Eq. (|34|). The specific heat calculated from Eq. ( |59"D and the self- 
consistent solution of the Bruekner equations is shown in Figure [X0| as a function of T at 
g = g c and g = 3.0. We have also calculated the specific heat by fitting a high-order poly- 
nomial to the internal energy obtained using quantum Monte Carlo simulations. As seen in 
Fig. [TU], the agreement between the analytical and numerical results is excellent for T < 0.5 J. 
The agreement is perfect at the critical point, while at g = 3.0 it is not so good. The reason 
for this can be traced to the fact that the value of triplet gap, A ~ 0.78 J, at T — 0, g — 3.0 
in our analytic method is slightly different from the numerical result Ao ~ 0.75J at zero 
temperature. This difference leads to a slightly different exponent for the specific heat. 

In the limit T <C J, g — > g c , the expression for the specific heat can be found explicitly 
by combining Eqs. (|l]), ((HD and (|59|): 

Cv = » A, (60) 



2tt7 2 v T 
where the universal function \1/ is given by 

1 r 00 dx\X\e Xl 



^>{x) 



2C(3) Jy( x) (e*i - If 



x 1 — y (x) + xy[x)- 



dx 



(61) 



The coefficient of \l/(x) in (^) has been chosen to be the specific heat of a single gap less 
Bose degree of freedom with dispersion u = 'jk in two dimensions. The value of ^(x) at 
T — >• is thus a measure of the effective number of such modes in the ground state. In 
the renormalized classical limit, it has the value \&(— oo) = 3, while at the critical point 
\P(0) = |\I/(— oo). In the quantum disordered limit, x — > oo, we find exponentially small 
specific heat: 

Cy w - — -re , x — > oo. 

Thus \I/(oo) = 0, corresponding to the absence of any gapless degrees of freedom. The 
limiting behavior of universal function \l/(x) is exactly the same as obtained one for the 
nonlinear cx-model (N — ► oo).0 

The fact that \E r (— oo) = 3 implies that our approach does not describe well the physics 
for g < g c , where the spin dynamics is well described by rotationally averaged spin-wave 
(Goldstone) fluctuations about a Neel-ordered ground state. There are only two gapless 
modes in the renormalized classical regime, while in our approach there are three, and the 
critical and Goldstone fluctuations are indistinguishable. The application of the Brueckner 
theory in this phase requires modifications; a nonzero sublattice magnetization, i.e., a single 
particle condensate at zero momentum, has to be introduced. 
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C. Uniform magnetic susceptibility 



Let us now calculate another important physical observable — the magnetic susceptibil- 
ity. If a uniform magnetic field of strength h is applied, then the term 



k 



(62) 



should be added to the Hamiltonian (|3p. 

Let us first to concentrate on the zero temperature properties. The second order correc- 
tion to the total energy of the system due to the perturbation (^) is given by the diagrams 
shown, in Figure |TT|(a), 



E {2) ~h 2 J2j^ [G n (K «)G„(k, w) - G? a (k, w)G«(k, to)] . 



(63) 



where G n (k,t) = -i(T(t koi (t)tl a (0))) and G a (k,t) = -i < T(il kQ (t)4 a (0)) > are normal 
and anomalous Green's functions at zero temperature which can be formally obtained by 
setting rik = in the expressions Q2"7D and (PB[). After a simple integration over u, one can 
show that = 0. Thus, the uniform magnetic susceptibility per spin, 



A 



ld 2 F 



2dh 2 



h=0' 



(64) 



is zero at T = 0. It is worth noting that the Bogoliubov transformation (0) does not change 



the operator of the magnetic moment M a = —ie a p^ X)k 



-le. 



k/37 Sk a k/3°k7- 



To calculate the uniform susceptibility at finite temperature, let us demonstrate that 
within the framework of the Brueckner approach the perturbation (|62|) splits a three de- 
generate triplet into excitations at U0k and loj, ± fi ■ h with magnetic moment fi = 1 at zero 
temperature. The correction to the self-energy due to the magnetic field h is given by the 



diagrams (b) and (c) in Figure ITT 



5SJ 7 (K) 



% hot ^Oi 



01 



r(K + Q) [G n (Q)G n (Q) - G a (Q)G a (Q)] 



d 3 Q 

(2tt) 3 



(65) 



- 5 / T 2 (K + Q)G n (K + Q - Qi)gn(Qi)Gn(Q) |^|^ ) 



:/3 7 



<9£ 



h=0 



duo 



w=0. 



ih a e a p.yZ] f 



-i 



The first diagram in Fig. |TT|(b) results in a trivial correction, ih a e a/ 3 y , the second and the 
third diagrams [the term linear in T(K) in eq.(|65|)1 give {i / 4)h a e a ^dT, / duo , and the diagrams 
in Fig. |TT](c) [the term quadratic in T(K) in the eq.(|65])] yield — (5/4)i/i a e Q ^ 7 <9£/<9u;. Further 
separation of a quasiparticle contribution in a Green's function for physical operators a^a 
with polarization a taken in spiral coordinates, 



uo — — E(k, to) + e ■ h ■ 



to — LUk + e ■ h' 
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with e = —1,0, 1, results in a simple Zeeman splitting of the elementary triplet excitation. 
Thus, instead of three degenerate triplets at there are three excitations at and uj^ ± h, 
i.e., for magnetic moment /x = 1. We can now calculate the magnetic susceptibility at 
finite temperature. Since the temperature correction to the renormalization constant is a 
very small, the main temperature correction to the self-energy is due to the second and 
third diagrams in Fig. 0(b). We use the conventional Matsubara technique (ui — > iu n ) to 
calculate this correction and obtain 



T,h 



-zAi a e a /3 7 ArpT2 2^ 



4TL 2 Y sinh 2 (cj q /2T)' 



(66) 



Using (]65|) and 
temperature: 



6|) one can find the quasiparticle magnetic moment renormalized due to 



Mp) 



z p r( P ,o) 



E 



4TL 2 ^ sinh i (cu q /2T) 



(67) 



Then the uniform magnetic susceptibility per spin can be easily obtained from 



I DM 



X 



2 dh 



h=0 



2 Bn 



(68) 



Figure [T2] shows the magnetic susceptibility as a function of temperature at g = g c and 
at g — 3.0, as obtained from Eq. (|Sg|) using the self-consistent solution of the Brueckner 
equations. In the same figure we also show our quantum Monte Carlo results. The agreement 
between the both methods is excellent at low temperature. The agreement at critical point 
is better than at g — 3.0 because the value of triplet gap, A ~ 0.78 J, at T = 0,g = 3.0 
in our analytic method is slightly different from numerical result Ao ~ 0.75. This difference 
leads to a slightly different exponent for the magnetic susceptibility. 

In the critical region, where T< J,^C g c , the main contribution to the integration in 
( [670 is due to small momenta, q ~ A/7, thus the susceptibility per spin can be written as 



T 



X 



27T7 : 



(69) 



with the universal function 



n(x) 



00 dx\X\t Xl 
y(x) (e^ - l) 2 



5 T 



,oc x\dxi 
•/(as) sinh 2 T 



(70) 



y(x)ctanh 



y{x) 



x 



1 — 25qT ( y(x)ctanh- x 



The constant 5 = Z C T C / (Sn^ 2 ) m 0.06 has been evaluated using the values of T c w 6.3 J, Z c w 
0.8,7 ~ 1.9 J taken from the self-consistent solution of the Brueckner equations at zero 
temperature, see Refs. ITRPOI. The universal function Q(x) has the limiting behavior 
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n(x)n^y (l-^x)(l- 25 TV5y ), x^O (71) 
n(x)w-|(l-^)(l + 25 A ), x^-oo (72) 
fi(x) ^xe- x (l + -), x-^oo, (73) 

which is exactly the same as predicted using the nonlinear a model (N — > 00)00 if correc- 
tions proportional to 5 are neglected. These corrections are subleading terms, but in the 
quantum critical and renormalized classical regimes they are parametrically larger then the 
subleading terms obtained for the a model (N — > oo). 



V. CONCLUSIONS AND DISCUSSION 

In conclusion, we have presented an effective, self-consistent diagrammatic approach to 
describe the properties of the quantum double-layer Heisenberg antifferrromagnet at low 
temperature and g > g c . To account for strong correlations between the elementary exci- 
tations (triplets), we applied the Brueckner approximation, treating the triplets as a dilute 
Bose gas with infinite on-site repulsion. We have also used a numerically exact quantum 
Monte Carlo simulation method to obtain non-perturbative results for comparison. We have 
calculated temperature dependent spectra (including the triplet gap), dynamic and static 
structure factors, the specific heat, and the uniform magnetic susceptibility. The agreement 
between the analytical and numerical results is excellent. 

Our analytical method involves approximations and an error of a few percent is always 
expected. The diagrammatic approach is valid for small on-site density of excitations because 
we neglected anomalous contributions in the vertex T(K). We emphasize, however, that the 
region of applicability of the present technique is surprisingly large. The reason for this can 
be traced to the fact that the density of triplets increases rather slowly as a function of 
temperature. We found that near the critical point: p T = 0.15 for T = 0.5J while p = 0.12 
at T = 0. It justifies the approach even for a relatively high temperatures, T ~ 0.5J, when 
9 > 9c- 

The Brueckner approach allows to calculate the elementary triplet spectrum at finite 
temperature, not only for small momentum but in the whole Brillouine zone. Comparing 
the analytical results for the dimensionless ratio -R(q) = S*(q)/Tx(q) with quantum Monte 
Carlo data, we demonstrated that the dynamic structure factor is peaked at u = u^, the 
broadening of the triplet excitation is negligible and the elementary triplet is well defined 
for g > g c ,T < 0.5 J. 

In the T J and g —>■ g c limit, our analytical results show that correlation length, the 
specific heat, and the magnetic susceptibility are universal functions of the zero tempera- 
ture gap, A , the spin-wave velocity, 7, and the temperature, T. The universal functions for 
the specific heat and the correlation length have the same limiting behavior as those pre- 
dicted using the nonlinear u-model (N — > 00), while our universal function of the magnetic 
susceptibility contains additional subleading terms. 

The agreement with the a model approach is due to the relative smallness of the quartic 
interaction, Eq. (^E|), for the two-layer Heisenberg model considered in this paper. In this 
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situation, the hard-core constraint is the most important and it is very natural that the 
results correspond to that of the a model. However, this situation is not general. There 
are many models where the quartic interaction is very important. It can even produce 
bound states of triplet spin-waves0 which effectively change the number of relevant degree 
of freedom. In this situation, one can expect a very substantial deviation from the simple 
nonlinear a-model behavior, while the above Brueckner approach can still be applied. An 
important example of such a system is the 2D J\ — J2 model where the singlet bound state 
has an extremely low energy.^ 

The Brueckner theory describes quite well the quantum critical and quantum disordered 
phases near the quantum critical point, as evidenced by the very good agreement with the 
quantum Monte Carlo results. However, the description of the renormalized classical regime 
is not good in this theory because the "Goldstone" regime with a Josephson length scale £j, 
where the spin dynamics is well described by rotationally averaged spin-wave fluctuations 
about a Neel-ordered ground state, is not present in our approach. In the Brueckner approach 
the critical and Goldstone fluctuations are indistinguishable. It is well known, however, that 
there are only two gapless modes in the renormalized classical regime, while in our approach 
there are three. The application of the Brueckner theory in this phase requires modifications; 
a nonzero sublattice magnetization N z , i.e., a single particle condensate at zero momentum, 
has to be introduced. We plan to consider this in a future study. 

The advantages of our analytic approach are that it is simple and captures the essential 
physics, as demonstrated here by comparing results with quantum Monte Carlo simulations 
as well as with predictions from the nonlinear a model. Obvious other applications of the 
Brueckner method include 2D Heisenberg models with frustration, Heisenberg ladders, as 
well one- dimensional spin chains at finite temperature. 

In addition to confirming the high accuracy of the Brueckner theory, the quantum Monte 
Carlo simulations performed here also show unambiguosuly that quantum critical behavior 
can be observed in the bilayer model at temperatures as high as T/J < 0.5. We also obtained 
an improved estimate of the critical coupling ratio; g c = 2.525 ± 0.002. 
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FIGURES 



FIG. 1. a) Equation for the scattering vertex T(K). b) Diagram for the self-energy S(K) 
corresponding to T(K). 

FIG. 2. Triplet gap, A/ J, as a function of coupling g = J±/J. Solid and dashed lines are the 
results of the self-consistent solution using the Brueckner approach and the Bruekner equations 
linearized in density, respectively. The dots with error bars are estimates obtained by dimer series 
expansions Jll 

FIG. 3. QMC results for the T = spin stiffness p s times the linear system size L, for 
L = 4, 6, . . . , 16, 20. The curves are quadratic fits to the QMC data (solid circles). The (neg- 
ative) slope of Lp s vs g increases with L. Statistical errors are typically of the order of the radius 
of the circles. 

FIG. 4. The uniform magnetic susceptibility, P er spin vs temperature, T/J, for different 
system sizes at a coupling g = 2.53. Statistical errors are much smaller than the sympols. 

FIG. 5. The uniform magnetic susceptibility, \J ■, P er spin vs temperature, T/J, for g = 2.52 
and 2.53. The lines are fits to the T/J < 0.17 data. 

FIG. 6. The ratio of the staggered structure factor and T times the staggered susceptibility 
for L = 32 (solid circles), L = 64 (open circles), and L = 128 (solid squares), for a coupling slightly 
below {g = 2.52) and above (g = 2.53) the critical coupling as a function of temperature, T/J. 
Statistical errors are at most of the order of the size of the symbols. The dashed lines indicate the 
constant value predicted within the nonlinear cr-model approach. 

FIG. 7. Deviation 5 A = (At — Ao)/J as a function of temperature, T/J, at the crititcal point, 
g = g c , and in the quantum disordered regime at g = 3.0. Solid lines are results of self-consistent 
solutions of the Brueckner equations linearized in density, dots are estimates extracted from Monte 
Carlo simulations for R(0) using Eq.(|57|), and dashed lines are predictions of the analytical Eq. ([48]) . 

FIG. 8. The dimensionless ratio R(0) = S(0)/[Tx(0)] as a function of temperature, T/J, for 
g = 2.52 and g = 2.53 obtained using quantum Monte Carlo simulation (dots). The solid line is the 
analytical result for g = g c obtained from the self-consistent solution of the Brueckner equations. 

FIG. 9. Elementary triplet spectrum, u>(q)/J, along the triangle in Brillouine zone, 
q = (0, 0) — (tt, it) — (it, 0) — (0, 0), at g = g c , T = 0.5 J. The solid line is the result of a self-consistent 
solution of the Brueckner equations, the dots are estimates extracted from Eq. ( |56|) and Monte 
Carlo results for R(q). 
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FIG. 10. Specific heat per spin, Cy, as a function of temperature T/J at critical point g = g c 
and g = 3.0. Solid lines are the results of self-consistent solution of Brueckner equations, dots are 
estimates obtained by Monte Carlo simulations. 

FIG. 11. a) Diagrams for the second order correction to the total energy due to a magnetic 
field perturbation. b),c) the diagrams for the correction to the normal self energy due magnetic 
field. 

FIG. 12. Magnetic susceptibility per spin, %«/, as a function of temperature, T/J, at g = g c 
and g = 3.0. The solid lines are results of self-consistent solutions of the Brueckner equations 
linearized in density. The dots are Monte Carlo results. 
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